NASA Technical Memorandum^. 10^3^ . 
ICOMP-93-37 



Multigrid Time- Accurate Integration of 
Navier-Stokes Equations 


Andrea Amone 

Institute for Computational Mechmics inJ*ropulsw^^ 
Lewis Research Center: ;zir ^ 

Clev ela nd, Ohio _ _ 

and University of Florence^ 

Florence, Italy 

and 


Meng-Sing Lion and I^uis A. Ppvinelli 

National Aeronautics and Space Administration 

Lewis Research Cent^ (nasa-tm-io6373) hultigrid N94-i7258 | 

Cleveland, Ohio — time-accurate integration of 

NAVIER-STOKES EQUATIONS (NASA) 

11 p Unclas - 

G3/34 0193181 


November 1993 


(W\SA 




MULTIGRID TIME-ACCURATE INTEGRATION OF NAVIER-STOKES EQUATIONS 

Andrea Amone* 

Institute for Computational Mechanics in Propulsion 
Lewis Research Center 
Cleveland, Ohio 44135 

and University of Florence 
Florence, Italy 

and 

Meng-Sing Liou^ and Louis A. Povinelli^ 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135 


Abstract 

EfBcient acceleration techniques typical of explicit 
steady-state solvers are extended to time-accurate 
calculations. Stability restrictions are greatly reduced by 
iwang of a fully implicit time discretization. A four- 
stage Runge-Kutta scheme with local time stepping, 
residual smoothing, and multigridding is used instead of 
traditional time-expensive factorizations. Some 
applications to natural and forced unsteady viscous flows 
show the capabili^ of the procedure. 


Introduction 

Recent progress in Computational Fluid Dynamics 
along with the evolution of computer performance is 
encouraging scientists to look at the details of flow 
physics more and more. There are a variety of practical 
applications where the unsteadiness of the problem can 
not be neglected (i.e. vortex shedding, natural 
unsteadiness, forced unsteadiness, aeroelasticity, 
turbomachinery rotor-stator interaction). Up to now, in 
several branches of engineering most of the analysis and 
Hpxi gning tools are based on a steady or quasi-steatfy 
assumption, even if the flow is known to be unsteady. 
Today, due to the improvement in computer resources, 
there is a strong interest in developing methodologies for 
efficient and reliable simulation of unsteady flow 
features. 

It is a common experience, while using time accurate 
explicit schemes, to be forced to choose the time step on 
the basis of stability restrictions. As a consequence, 
unless the problem is a very high frequency one, the 
number of time steps to be performed is much higher 
than the one required for time accuracy. 
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By tnftans of some implicit factorization, most of the 
stabili^ restrictions can be removed, but the work 
required at each time step grows rapidly with grid 
dimension and complexity of the flow equations. In 
addition, several boundary conditions can be difficult to 
treat in a fully implicit way. 

In viscous flow calculations, the grid is stretched 
close to the shear layer and the characteristic time step 
varies several orders of magmtude inside the 
computational domain. Even if in practical applications 
the characteristic time step of the core-flow region is 
comparable with the one suggested by accuracy, close to 
the walls the time step restrictions are very severe. 
Therefore highly vectorizable schemes with less stability 
restrictions on the allowable time step would be an 
interesting combination. 

Explicit schemes with accelerating techniques have 
proven to be very effective for solving steady 
problems‘’2"3. Unfortunately, the computational 
efficiency of those time-marching solvers is achieved by 
sacrificing the accuracy in time. In this paper, we present 
a procedure to show that the conventional steady-state 
acceleration techniques, specifically the multigrid 
tftr hni qiiKs, can Still be applied to unsteady Navier- 
Stokes problems as well, while still achieving efficiency. 
The basic idea is to reformulate the governing equations 
so that they can be handled by an explicit accelerated 
scheme'*. If the time discretization is made implicit, 
stability restrictions are removed and accelerating 
tp^hniqiip-fi can be used instead of traditional time- 
expensive factorizations {i.e. ADI, LU). 

As one of the final goals of the present research will 
be the study of unsteady phenomena in turbomachinery 
conqwnents, such as rotor-stator interaction and stage 
analysis, we implemented the techmque in the 
TRAF(2D/3D) codes^-^. These two- and three- 
dimensional solvers were developed during a joint 
project between the University of Florence and NASA 
Lewis and were designed for tuibomachinery blade row 
analysis. 



The procedure is validated by applying it to some 
examples of natural and forced unsteady two- 
dimensional viscous flows. 


Governing Equations 


Let t,p, u, V, p, T, E, and H denote respectively time, 
densi^, the absolute velocity components in the x and y 
Cartesian directions, pressure, temperature, specific tot^ 
energy, and specific total enthalpy. The two- 
dimensional, unsteady, Reynolds-averaged Navier-Stokes 
equations can be written for a moving grid in 
conservative form in a curvilinear coordinate system ij 
as. 


Xy,-2flVy+^Ul-i- V,.) 

Xxy — Tyx “^Vjt) (7) 

Py=UX„+VXyy+kT, 
P,=UXy.+VXyy+kTy 

and the Cartesian derivatives of (7) are expressed in 
terms of ^-, and 7 - derivatives using the chain rule, i.e., 

«x=i,u( + Tj^un (8) 

The pressure is obtained from the equation of state. 


dt dr] dr] 


where. 





+^yP 

pHU + ^,p 


p=pRT (9) 

According to the Stokes hypothesis, 2. is taken to be 
-2p/3 and a power law is used to determine the 
molecular coefficient of viscosity p as function of 
temperature. The edtly-viscosity hypothesis is used to 
account for the ^ect of turbulence. The molecular 
viscosity p and the molecular thermal conductivity k are 
replaced with, 

p=p,+p, ( 10 ) 


The contravariant velocity components of eqs. ( 2 ) are 
written as, 

u =4,+^,u+4yV, V = rj, + T],u + Tj,v (3) 

and the transformation metrics are defined by, 

4y=Jxn, 4,=-x,^,-y,4y 
x]y=Jx(, n,=-x,Ti,-y,T}^ (4) 

where the Jacobian of the transformation J is, 

J~'=xiy^-x„y^ (5) 

The viscous flux terms are assembled in the form, 

0 ] [ ^ 

Tjy ^xTxx~^ Vy Txy 

4i Xyx 4y Xyy 7t Xyx Vy Xyy 

4A+4yPy J [nJ, + lyPy 

where, 

r*r =2pux +2.{ux+v^ 





where cp is the specific heat at constant pressure, Pr is 
the Prandtl number, and the subscripts / and t refer to 
laminar and turbulent quantities respectively. The 
turbulent quantities and Pr, are computed using the 
two-layer mixing length model of Baldwin and Lomic^. 


aatial Discretization and Artificial Pissipation 


Traditionally, using a finite-volume approach, the 
governing equations are discretized in space starting 
from an integral formulation and without any 
intermediate mapping. In the present work, due to the 
large use of eigenvalues and curvilinear quantities, we 
found it more convenient to map the Cartesian space 
(x,y) in a generalized curvilinear one (^, 1 ]). In the 
curvilinear system, the equation of motion ( 1 ) can be 
easily rewritten in integral form by means of Green's 
theorem and the metric terms are handled following the 
st^dard finit^volume formulation. A cell-centered 
scheme is used to store the flow wi^les. On each cell 
face the convective and difliisive fluxes are calculated 
after computing the necessary flow quantities at the face 
center. Those quantities are obtained by a simple 
averaging of adjacent cell-center values of the dependent 
variables. 
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In viscous calculations, dissipating properties are 
present due to difiusive terms. Away from the shear 
layer regions, the physical difiiision is generally not 
sufficient to prevent the odd-even point decoupling of 
centered schemes. Thus, to maintain stability and to 
prevent oscillations near shocks or stagnation points, 
artificial dissipation terms are also included in the 
viscous calculations. Equation (1) is written in semi- 
discrete form as, 

^+C(Q)-D(Q)=0 (14) 

where the discrete operator C accounts for the physical 
convective and difiusive terms, while D is the operator 
for the artificial dissipation. The artificial dissipation 
model used in this paper is basically the one originally 
introduced by Jameson, Schmidt, and TurkeF. In order 
to minimize the amount of artificial difiiision inside the 
shear layer, the eigenvalues scalings of Martinelli and 
Jameson^, and Swanson and Turicel^ have been used to 
weight these terms. The quantity D(Q) in eq. (14) is 
defined as. 


D(Q) -D^f+iy,-Dl)Q (15) 

where, for example, in the # curvilinear coordinates we 
have, 

DfQ = Vi(A UI/2J 

DfQ = VfAfQij (16) 


ij are indices associated with the ^,tj directions and 
p^, 2 I 4 are forward and backward difference operators 
in the ^ direction. The variable scaling factor A is 
defined as. 


and. 


A^ — 


(17) 

(18) 


<pf=l + 



(19) 


where and A^ are the scaled spectral radii of the flux 
Jacobian matrices for the convective terms. 


( 20 ) 


and a is the speed of sound. Note that the effect of the 
grid motion is accounted for in ( 20 ) through the 
definition of the contravaiiant components of velocities 
of (3). The exponent tr is generally defined by 0<a<l , 
and for two-dimensional applications, a value of 2/3 
gives satisfactory results. The coefficients and 
use the pressure as a sensor for shocks and stagnation 
points, and are defined as follows. 


=K!^^MAX( ( 21 ) 


Pi-ij ^ Pij^Pi^tj 

Pi-ij 




( 22 ) 

(23) 


where typical values for the constants and are 
1/2 and 1/64 reqtectively. For the other direction, rj, the 
contribution of dissipation is defined in a similar way. 
The computation of the dissipating terms is carried out 
in each coordinate direction as the difference between 
first and third difference curators. Those operators are 
set to zero on solid walls in order to reduce the global 
error on the conservation property and to prevent the 
presence of undamped modes®-*®. 


Boundary Conditions 

In cascade-like configurations there are four different 
types of boundaries; inlet, outlet, solid wall, and 
periodicity. According to the theory of characteristics, 
the flow angle, total pressure, total temperature, and 
isentropic relations are used at the subsonic-axial inlet, 
while Ae outgoing Riemann invariant is taken from the 
interior. At the subsonic-axial outlet, the average value 
of the static pressure is prescribed and the density and 
components of velocity are extrapolated. 

On the solid walls, the pressure is extrapolated from 
the interior points, and the no-slip condition and the 
temperature condition are used to compute density and 
total energy. For the calculations presented in this 
paper, all the walls have been assumed to be at a 
constant temperature equal to the total inlet one. 

Cell-centered schemes are generally implemented 
using phantom cells to handle the boundaries. The 
periodicity is, therefore, easily overimposed by setting 
periodic phantom cell values. On the iMundaries where 
the grid is not periodic, the phantom cells overlap the 
real ones. Linear interpolations are then used to 
compute the value of the dependent variables in phantom 
cells. 
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Basic Time-Stepping Scheme and Acceleration 
Techniques for the Steady Problem 

The system of the differential equation (14) is 
advanced in time using an explicit four-stage Runge- 
Kutta scheme until the steady-state solution is reached. 
A hybrid scheme is implemented, where, for economy, 
the viscous terms are evduated only at the first stage and 
then fiozen for the remaining stages. If / is the index 
associated with time we wiU write it in the form. 


dr = 






1 


where yis the specific heat ratio and. 


(27) 


(28) 


e'"=e' 


(2«) 

K, being a constant whose value has been set equal to 2.5 
based on numerical experiments. 


(24) 

111 , 
ai=—, 02=—, 03=—, 04 =i 

4 i 2 


Residual Smoothing 

An implicit smoothing of residuals is used to extend 
the stability limit and the robustness of the basic scheme. 
This technique was first introduced by Lerat** in 1979 in 
conjunction with Lax-Wendroff ^pe schemes. Later, in 
1983, Jameson^ implemented it on the Runge-Kutta 
stepping scheme. In two dimensions, the residual 
smoothing is carried out in the form. 


where the residual R (Q) is defined by, 

R(Q)=AtJ[C(Q)-D(Q)] (25) 

Good, high-frequency damping properties, important for 
the multigrid process, have been d>tained by perfomiing 
two evaluations of the artificial dissipating terms, at the 
first and second stages. 

In order to reduce the computational cost, three 
techniques are employed to speed up convergence to the 
stea<fy state-solution. These techniques: 1) local time- 
stepping; 2) residual smoothing; 3) multigrid; are briefly 
described in the following. 


{l-fi,V'fAi){j-fi,V,A,)R=R (30) 

where the residual R includes the contribution of the 

variable time step and is defined by (25) and R is the 
residual after a sequence of smoothing in the and 7 , 
directions with coefficients , and . For viscous 
calculations on highly stretched meshes the variable 
coefficient formulations of Martinelii and Jameson^ and 
Swanson and Tuikel^ have proven to be robust and 
reliable. In the present paper, the expression for the 
variable coefficients P of (30) has been implemented as 
follows. 


Local Time-Stepping 

For steady state calculations, a faster expulsion of 
disturbances can be achieved by locally using the 
maximum allowable time step. In the present work the 
local time step limit At is computed accounting for both 
the convective (At^ and diCRisive (At^ contributions as 
follows. 


f 

At =co 

V 




(26) 


where is a constant usually taken to be the Courant- 
Friedrichs-Lewy (CFL) number. Specifically, for the 
inviscid and viscous time step we used. 


P^=MAX< 



CFL Jlj 
CFL* 





P,=max\ 



CFL Xr, 
CFL* X^+Xr, 


V 



J 



(31) 


where the coefficients 0 ^ and <P,^are the ones defined in 
eqs. (19), and CFL, and CFL* are the Courant numbers 
of the smoothed and unsmoothed scheme respectively. 
For the hybrid four-stage scheme we used CFL=5, and 
CFL*=2.5. 
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Multigrid 

This technique was developed in the beginning of the 
1970s for the solution of elliptic problems^^ and later 
was extended to time-dependent formulations*’^. The 
basic idea is to introduce a sequence of coarser grids and 
to use them to speed up the propagation of the fine grid 
corrections, resulting in a faster expulsion of 
disturbances. In this work, the Full Approximation 
Storage (FAS) schemes of Brandt*^ and Jameson* is 
used. 

Coarser, auxiliary meshes are obtained by doubling 
the mesh spacing and the solution is define! on them 
using a rule which conserves mass, momentum, and 
energy, 

=z(j-^q) (32) 

2h h 

where the subscripts refer to the grid spacing, and the 
sum is over the eight cells which compose the 2h grid 
cell. Note that this definition coincides with the one 
used by Jameson* when the reciprocal of the Jacobians 
are replaced with the cell volumes. To respect the fine 
grid approximation, forcing functions P are defined on 
the coarser grids and added to the governing equations. 
So, after the initialization of Q 2 ), using eq.(32), forcing 
functions P 2 h are defined as, 

P2H=lRdQ,)-R2H(^2H) ( 33 ) 

and added to the residuals R 21 , to obtain the value Pi 2 h 
which is then used for the stepping scheme. 

Rsh = Rlh( Q21J +P2h (34) 

This procedure is repeated on a succession of coarser 
grids and the corrections computed on each coarse grid 
are transferred back to the finer one by bilinear 
interpolations. 

A V-type cycle with subiterations is used as a 
multigrid strategy. The process is advanced from the 
fine grid to the coarser one without any intermediate 
interpolation, and when the coarser grid is reached, 
corrections are passed back. One Runge-Kutta step is 
performed on the h grid, two on the 2h grid, and three on 
all the coarser grids. 

For viscous flows with very low Reynolds number or 
strong separation, it is important to compute the viscous 
terms on the coarse grids, too. The turbulent viscosity is 
evaluated only on the finest grid level and then 
interpolated on coarse grids. 

On each grid, the boundary conditions are treated in 
the same way and updated at every Runge-Kutta stage. 
For economy, the artificial dissipation model is replaced 
on the coarse grids with constant coefficient second- 
order differences. On coarse grids, the turbulent 


viscosity is evaluated by averaging the surrounding fine 
grid values. 

Reformulation of the Governing Equations 

Explicit Runge-Kutta schemes in conjunction with 
residud smoothing and multigrid have proven to be very 
efficient for steady problems, however those time- 
dependent methods are no longer time accurate. As 
shown by Jameson'* for the Euler equations, the system 
of (1) can be reformulated to be handled by a time- 
marching steady-state solver. The equations (1) and (14) 
are rewritten in a compact form as, 

^=-mQ) ( 35 ) 

at 

where is the residual which includes convective, 
diffiisive, . and artificial dissipation fluxes. By the 
introduction of a fictitious time r the unsteady governing 
equations can be reformulated and a new residual iV* 
defined as, 

(36) 

dr dt 

now r is a fictitious time and all the accelerating 
techniques developed in stea<fy-state experiences, can be 
used to efficiently reduce the new residual iH* , while 
marching in x. Following the approach of Jameson'*, 
derivatives with respect to the real time t are discretized 
using a three-point teckward formula which results in an 
implicit scheme which is second order accurate in time, 

I + ^(Q^‘)=pr* (or‘) (37) 

dr 2 At 

where the superscript n is associated with the real time. 
Between each time step the solution is advanced in a 
non-physical time r and acceleration strategies like local 
time stepping, implicit residual smoothing, and 
multigridding are used to speed up the residual SR* to 
zero to satisfy the time-accurate equations. 

The time discretization of (37) is fully implicit, 
however, when solved by marching in z, stability 
problems can occur when the stepping in the fictitious 
time rexceeds the physical one. This generally occurs in 
viscous calculations where core-flow cells are much 
bigger then those close to shear-layer. Based on a linear 
stability analysis of the four-stage scheme of (24) applied 
to (37), the stepping in x must be less then 2/3 CFL'At. 
The time step zlrcan then be corrected as follows, 
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( 




Ax=MlN\ 


Ax, 




At 

yt-l ^ CFL 
2 CFL* > 


(38) 


where the contribution of tt^ muitigrid speed-up is 
included through 2***-^, m being the total numto of grids 
used in the multigrid process. After limiting the time 
step Ax with (38) the scheme becomes stable and the 
physical time step At can be chosen safely only on the 
basis of the accuracy requirement 

At the end of each time step in real time, the time 
derivative ^Q/dt is updated and a new sequence in the 
fictitious time x is started. From 10 to 20 multigrid 
cycles are typically needed between time steps. 

To provide a good initialisation of the solution at the 
new time step, a three-point backward formula is used as 
a predictor. 


Q- (39) 

where Q* is the predicted value of 0^^, 

We stress that, using scheme (37), the modifications 
to turn an existing steady-solution solver to a time 
accurate one are quite simple. The time derivative dQ/dt 
can be introduced as a source term to be included in the 
new residual X*, and the time step is corrected using 
eq. (38) to make the scheme stable. 


Results and Discnsslons 


In this very first part of the research project the 
procedure is validated in two-dimensions. Three test 
cases are presented. Firstly, a vortex shedding over a 
row of circular cylinders in a laminar regime is 
examined. The interest being mostly in the flow 
periodicity and in the prediction of the Strouhal number. 
As a second application of natural unsteadiness, a shock 
buffeting over a row of bicircular airfoils will be 
discussed. Finally, the last application is related to 
forced unsteadiness in turbomachines and simulates the 
effect of passing stator wakes on a rotor blade. 

Row of Circular Cylinders 

This test is intended to predict the natural vortex 
shedding past a cylinder, A row of circular cylinders in 
a laminar regime is studied for an inlet flow condition of 
Mach number 0.2 and Reynolds number of 1000. 
Calculations were performed on a 257x49 elliptic C-type 
grid. The distance between the cylinders is five times the 
cylinder diameter. 

Figures 1(a) and (b) report the evolution in time of 
the flow angle and velocity components (phase plot) at a 


point in the wake close to the cylinder. The time history 
refers to four cycles of oscillations after a periodic flow 
condition is reached. The very pericxiic betaviour of the 
flow is evident and proves the rc^ustness and accuracy of 
the scheme. The time stq> for those calculations was set 
to have 40 divisions over a cycle. This corresponds to a 
local Courant number between three (far field) and four 
hundred (boundary layer). 

The computed Strouhal number based on the inlet 
velocity is about 0.2 and agrees well with the 
experimental value of 0.21. 

Figure 2 reports the instantaneous particle traces in 
nine instants over a cycle (the tenth position would be 
eqxiivalent to the first). The shedding of the vortex is 
very evident as well as the mechanism of their formation 
with a vortex merging between instants 1 and 2, and 5 
and 6. 



a) flow angle evolution 



b) velocity component evolution 


Fig 1. Unsteady flow past a circular cylinder 
(M=0.2,Re=10a). 
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Time=l Time=2 




Time=5 Tiine=6 



Time=7 Time=8 



Timc=9 



Fig. 2. Instantan eous paiticle traces for the circular 
cylinder row test case 


Far dW 3 y from the nominal condition, end^vall effects 
become important especially at high Mach numbers and 
so the comparison with a row of airfoils becomes not too 
meaningful. 



a) pressure coefficient distribution 



Shock Buffeting Over a Bicimilar Airfoil 

S tarting from about 1976 several experiments and 
calctilations were carried out on shock buffeting over a 
bicitcular-arc airfoil. The experiments*^'*'* were carried 
out at NASA Ames in a wind tuimel designed for this 
purpose. At a Reynolds number of 7 million, experiments 
suggest buffeting at a free stream Mach number in the 
range from 0.74 to 0.76. In agreement, the present 
caloilations indicate natural unsteady flow at Mach 0. 75 
On the contraiy, while, the flow is experimented to be 
unstearfy up to a Mach munber of ft 78, the calculation 
still shows some unsteadiness up to a free stream Mach 
number of 0.83. The nominal test Mach number for the 
experiments is ft 775 and the wind tuimel endwalls were 
designed to minimize wall effects for this flow condition. 



b) instantaneous Mach number contours 

Fig. 3. Shock buffeting over a row of bicircular-arc 
airfoils (A/=ft775, Re=7xl0^) 


7 


The reduced frequency of the experiment is roughly 
0.5. Steger**, with an isolated airfoil predicted about 
0.41, while the TRAF2D code suggested 0.42 for an 
airfoil distance of ten times the axial cord. If the airfoils 
are clustered to a distance of four times the axial cord, 
the reduced frequency rises to 0.47, once again 
suggesting some influence of existing walls on the 
buffeting fiequent^. 

Instantaneous Mach number contours ate reported in 
Fig. 3 along with the range of pressure distrilnition. For 
that high Rqmolds number the Courant number based on 
forty^ divisions within a cycle was between one (far field) 
and three thousand (boundary layer). The grid used is an 
H-^pe (153x97) and a buffeting cycle requires about 8 
minutes on the NASA Lewis Cray Y-MP. 



Passing Wakes Effects on a Rotor Blade 

As a preliminary application to unsteady effects in 
turbomachinery, a rotor configuration with incoming 
moving wakes was studied. A H-type grid (see Fig. 4(a)) 
was selected in order to minimize the incoming wake 
smearing due to grid coarsening. The wake is simulated 
with a loss in total pressure and an alteration in the 
velocity directiotL Figure 4(b) shows the isentropic 
Mach number range on the blade surface. The wind 
tunnel data^^ without wake effects are also reported. The 
instantaneous Mach number contours are given in Fig. 4 
(c) in terms of four instants over a qele. 



b) isentropic Mach number distribution on the blade 
surface 


Ttmc=l 



Timc=2 



Timc==3 


Timc=4 



c) instantaneous Mach number contours 


Fig. 4. Passing wakes effect on a rotor blade (M2,^=.S7, Re 2 = 8xl(fi) 
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Conclusions 


The use of explicit accelerated schemes has been 
extended to time-accurate Navier-Stokes calculations. 
Particularly the use of efficient and highly vectoiizable 
techniques such as multigridding is proposed in 
conjunction with a hilly implicit time discretization. A 
preliminary validation indicates that the approach is 
robust and efficient. According to the propos^ method, 
the modifications to be made on a time- marching 
accelerated solver in order to make it time accurate are 
very simple. 
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